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A SURVEY OF THE METHODS DEVELOPED 
FOR THE INVERSION OF THE RADIATIVE 
TRANSFER PROBLEM FOR PLANETARY ATMOSPHERES 


by 

Barney J. Conrath 
Goddard Space Flight Center 


SUMMARY 

Five of the methods developed for obtaining the thermal profile of a 
planetary atmosphere from remote radiometric measurements are discussed 
and compared. An extension of the "variable slab" method of King is used in 
making an inversion calculation for a model atmosphere. The question of the 
stability of the solutions of the basic integral equation of the problem and its 
implications on the accuracy of the derived temperature profile are con- 
sidered. It appears that some smoothing will have to be introduced, regard- 
less of the method used, in order to achieve a stable solution from a given set 
of observational data containing random errors. A considerable amount of 
work remains to be done, including extension of the variable slab method to 
arbitrary band transmissions and the performance of detailed error analyses 
of the various inversion methods. It also seems desirable to examine the 
problem from the point of view of the uses to be made of the derived thermal 
structure in order to ascertain whether the smoothed solutions retain suffi- 
cient physical information to be of value. 
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A SURVEY OF THE METHODS DEVELOPED 
FOR THE INVERSION OF THE RADIATIVE 
TRANSFER PROBLEM FOR PLANETARY ATMOSPHERES 


by 

Barney J. Conrath 
Goddard Space Flight Center 


INTRODUCTION 

The possibility of extracting the temperature profile of a planetary atmosphere from obser- 
vations of the emerging specific intensity, either as a function of viewing angle (Reference 1) or 
as a function of frequency (Reference 2), has been recognized for a number of years. In this 
survey the various methods which have been devised for obtaining solutions to this so-called 
’'inversion problem” will be summarized. When possible the methods will be compared, and the 
strengths and weaknesses of each method will be examined. The question of stability of solutions 
will be discussed, and an attempt will be made to define the areas of work remaining to be done 
on the problem. 


STATEMENT OF THE PROBLEM 

If the thermally emitted radiation emerging from the top of a planetary atmosphere is ob- 
served in a spectral region of relatively strong atmospheric absorption, it is evident that very 
little of the radiation originating from deep within the atmosphere will contribute directly since 
that part of the atmosphere simply cannot be "seen”. The regions near the top of the atmosphere 
cannot be seen due to their transparency. Hence, there must be a region of maximum contribu- 
tion to the outgoing radiation located some place between the two extremes. From this line of 
reasoning it is evident that we should be able to learn something about the temperature of the at- 
mosphere in these regions of maximum contribution by observing the emerging radiation in an 
absorption band. This method has recently been applied to the 15 (jl C0 2 absorption band in the 
earth’s atmosphere, using data from the Tiros VII meteorological satellite (Reference 3). By 
extending the observations to cover several bands of different atmospheric absorption it should 
be possible to obtain information on the temperature in several altitude regions of the atmosphere. 
Similar results should be obtainable by looking at several viewing angles, and hence through dif- 
ferent absorber masses, using only one spectral band. 
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To put these physically intuitive concepts on a quantitative basis, it is necessary to examine 
the solution of the equation of radiative transfer. For an atmosphere with sufficient optical thick- 
ness such that essentially no radiative contribution from the surface of the planet is received, the 
specific intensity emerging from the top of the atmosphere, in zenith direction 9 = cos' 1 /z, in a 
narrow frequency band centered at wave number v, can be written in the plane parallel atmosphere 
approximation as 


1(0, /X, v) 


dt du 

B ( u ) du 1 T 


(la) 


where t is the band transmission from the level of normal absorber concentration u to the top of 
the atmosphere, and B(u) is the Planck function averaged over the band interval. The normal ab- 
sorber mass is defined by the relation 


u 



(lb) 


where p is the density of the optically active gas and z is the altitude. 

It is sometimes convenient to write the integral in Equation la in terms of altitude or logarith- 
mic pressure scaling, i.e., for normal viewing: 


/•CD /• CC y%03 

B ( u ) fu du = B ( z ) fi dz = " B (P) 5(logp) d ( lo gP) 

Jq Jq Jo 


(lc) 


Thus, for an optically thick atmosphere, dt/dz or dt/d(logp) can be regarded as an normalized 
"weighting function" which gives essentially the relative contributions of each layer of the atmos- 
phere to the emerging intensity. Calculation of such weighting functions for regions of relatively 
high atmospheric absorption show more or less well defined maxima, the altitudes of which de- 
pend on the opaqueness of the spectral regions considered (cf. calculations of dt/d(logp) vs. p by 
Wark and Yamamoto for the 15 m C0 2 regions for conditions applicable to Mars in Reference 4). 
Thus, the quantitative formulation agrees with the physically intuitive line of reasoning given above. 

Let us assume that we have available measurements of 1(0, v) as a function of wave number. 
Assume for simplicity that the observations are made looking straight down (m = 1) . Then letting 
-<9t(u, v)/du s K(u, v) we have from Equation la, 


1(0, v) 


J 0 


K(u, v ) B(u) du 


( 2 ) 
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We wish to extract the temperature profile which is equivalent to finding the Planck function B( u ). 
The problem is mathematically one of solving a Fredholm integral equation of the first kind, of 
which Equation 2 is an example. The measured quantity 1(0, v) can be regarded as an integral 
transform of the Planck function with K(u, v) as the transform kernel. 

For the special case of strictly monochromatic radiation for which the mass absorption coef- 
ficient k does not vary along the path length, the intensity in direction 6 - cos " 1 m in the plane- 
parallel atmosphere approximation becomes 

l(o. Jr) = j B(U) e-"-"* J du . (3) 

In this case 1(0, k//x) is simply the Laplace transform of B and the solution can be formally written 


B(u) 



(4) 


where the path of integration is parallel to the imaginary axis and to the right of all singularities 
in the integrand. It should be noted that determination of B from measurements of I vs. v and 
from I vs, n are formally equivalent. 

Despite the formal simplicity of the mathematics of the problem, the solution of Equation 2 for 
B to a reasonable degree of accuracy for a given discrete set of measurements of X, containing ex- 
perimental errors, turns out to be a formidable task. Part of the reason for the difficulty is as- 
sociated with the fact that an inversion of Equation 2 in effect amounts to performing a differentia- 
tion of the observational data. We will return to this aspect of the problem in a later section after 
first considering in the next section several methods of solution which have been proposed. 


METHODS OF SOLUTION 
Perturbation Method 

In this method, proposed by Kaplan (Reference 5), a number of representative model atmos- 
pheres are used to calculate the emerging intensity in several narrow spectral bands located in a 
region of strong absorption such as the 15 m C0 2 region. Simultaneous intensity measurements 
made in the same spectral bands are used to calculate corrections to a ’’best guess” model atmos- 
phere chosen from the representative set on the basis of season and location. The calculation 
scheme used is 


1(0. - Io (o. V,) 

lo (0. "ij 




(5) 
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where I 0 is the intensity calculated from the model atmosphere, l(vj is the observed inten- 
sity and the AT ’s are the temperature corrections. One disadvantage of this system is the neces- 
sity of dividing the atmosphere up into arbitrarily chosen layers, the choice of which appears to 
have a considerable effect on the accuracy of the solution. An error analysis of this method by 
Drayson (Reference 6) indicates that the solutions are relatively insensitive to systematic errors 
in the measured intensity, but non- systematic errors of the order of 1 % result in totally unreal- 
istic values for the At t s of the order of thousands of degrees, indicating a lack of stability in the 
solutions. 


Analytic Method 

This method is based on the assumption that the Planck function can be well represented by a 
truncated expansion of the form 


B(u) 


n 



i ^1 


(6) 


where the F r s are a known set of functions. If the F’s are a complete set of orthogonal functions, 
then as n approaches infinity Equation 6 becomes exact. When Equation 6 is substituted into Equa- 
tion 2, with intensity observations atn different frequencies, we obtain the following set of equations: 


1 ( 0 , -) 



u ) Fj (u) du ; 


i - 1 , 2, * * * n . 


( 7 ) 


Since the integrals involve only known functions, they can be evaluated, yielding a linear system 
of n equations for the n a/s. 

This method has been applied by Yamamoto (Reference 7) to several different model atmos- 
pheres. Calculations were made with three different function sets: power series, Legendre 
polynomials, and Chebyshev polynomials. The calculations indicate that the function class giving 
the best representation differs from one model to the next. Yamamoto suggests that a statistical 
study, using many atmospheric soundings, might be made in order to determine the best function 
class to use for a given latitude and season. 

One obvious advantage of this method is the fact that it is not necessary to divide the atmos- 
phere into arbitrary layers. The chief disadvantage seems to be the problem of choosing the best 
set of functions for approximating the Planck function, since the expansion will always be limited 
to a small number of terms because of the limited number of observations available. 
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Imbedded Source Method 


This is a method proposed and developed by King (References 8 and 9). A grey plane- parallel 
atmosphere with intensity observations as a function of viewing angle was assumed, permitting the 
use of Equations 3 and 4. It was recognized that Equation 4 could be used if a reliable interpola- 
tion formula for the discrete observations 1(0, /x) could be found. 

In an effort to arrive at an optimum interpolation formula, the common physical factor in- 
fluencing both B(u) and 1(0, /x), namely the absorption of radiation externally incident on the at- 
mosphere, was considered. In using this approach, the integral equation form of the radiative 
transfer equation, 


dl ( t, m) 
M &? 


I (r,/i) 




I(t,/x')c^x # + 


dF 

dr 


was solved, assuming a source function of the form 


dF 

dr 


S(r ) 


F(0) 

M 0 6 


( 8 ) 


(9) 


where r is the normal optical depth in the infrared r = k ir u. Physically, this is a source due to 
a flux ttF(O) incident normally downward on the top of the atmosphere, attenuated exponentially as 
it penetrates the atmosphere. In this case /x 0 is equivalent to the ratio k ir /k visible , due to the 
definition of r as the infrared optical depth, and should not be confused with the usual viewing 
angle definition associated with /x 0 . Equation 8 has been solved with source function (Equation 9) 
by making the identification of the formal mathematical equivalence of the resulting integr ©differ- 
ential equation with that for conservative isotropic scattering (References 10 and 11). The result- 
ing solution for the intensity at the top of the atmosphere as a function of viewing angle is 


F( 0) h (m) h (^o) 

1(0, fi) - 4 ^ + Mo 


( 10 ) 


where H(/x) is the Chandrasekhar H- function defined as the solution of a certain nonlinear integral 
equation (Reference 11). The H-functions can be expressed to an approximation as quotient poly- 
nomials of the form 


H(aO 



n k h 

i = i 

IE l 1 + k a m) 

a= 1 


( 11 ) 


where the k T s are roots of an associated indicial equation. 
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This led to the idea that 1(0, /x) could best be represented by a quotient polynomial for inter- 
polation purposes of the form 


1(0, aO 


a n^ n + a n-l M * 1-1 + ' ' • 
b n + b n-l M "" 1 + • * • 



( 12 ) 


When attempts were made to fit this form to realistic data, it was found that singularities were 
obtained between the data points. This, of course, was physically inadmissible since 1(0, ^ ) is 
bounded. 

Since the source function assumed above was inherently smooth, the opposite extreme 
s ( t ) - $ (r - was tried next. This resulted in an interpolation formula which was of quasi- 
quotient polynomial form with the numerators given by exponentials. Apparently the method has 
not been carried much beyond this point. 

The advantages and disadvantages of this method are hard to ascertain since it has not been 
fully developed and extended to band transmissions. 


Fixed Slab Method 

In its simplest form this method replaces the integral in Equation 2 with a finite sum con- 
taining a number of terms not exceeding the number of data points available for 1(0, ^). In this 
way a system of linear equations is obtained in the form 


I (0, v i ) - y K (uj , B (u Auj ; i - 1,2, • • • , m; m > n . (13) 

j -l 

When m = n , a linear system of n equations in the n unknowns b(uj ) is obtained. If the number of 
observations exceeds the number of terms in Equation 13 (m > n) , then a least squares solution can 
be obtained. This method has been used by Wark (Reference 12), modified such that the quantity 
solved for is dB/dT . While the above description has assumed observations of intensity as a func- 
tion of frequency, the same type of formulation can be used for intensity as a function of viewing 
angle. 

A disadvantage of this method is the necessity of arbitrarily dividing the atmosphere into 
slabs in order to specify the u . T s in Equation 13. 


Variable Slab Method 

The variable slab method is the most recently developed of the methods discussed in this 
paper (Reference 13). A plane- parallel atmosphere with a simple exponential transmission is 
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assumed so that Equation 3 can be used. The method can be applied to either the I vs. v or the 
I vs.m type of observation, but for simplicity we shall consider only the former case in detail. 

Integrating Equation 3 by parts gives 


1(0, v) - B( 0) 



e - k <v)u du ^ 


(14) 


where B(0) is the Planck function at the top of the atmosphere. Next, the atmosphere is broken up 
into slabs of constant b, giving 


dB(u) 

du 


2 (AS). S (u - Uj) 
j =1 


(15) 


Substitution into Equation 14 yields 


1(0, V) - B(0) = 22 (AB) j e" k(v)Ui . (16) 

J =1 

Now instead of arbitrarily assigning the values of u . , 2n observations are used 
to find the n values of (AB) j and the n u^s. The 2n values of k(v) are required to obey the relation 


kU) 

= 1 + 1; l = 0, 1, 2, ••• , 2n - 1 , (17) 


where v 0 is an arbitrarily chosen base frequency. Introducing the following definitions, 


I (o , Vj) - B(0) 


(AB), e 


-k(^ 0 )u . 


(18) 


Equation 16 becomes 


( v o)u . 




l = 0, 1, 2, * •• , 2n - 1 . 


j =1 


(19) 
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Thus, the problem becomes one of solving 2n nonlinear equations for the n b^s and the n x/s. An 
algorithm for solving a set of equations of this type is given by Chandrasekhar (Reference 11) in 
connection with a discussion of quadrature formula. A set of n numbers C k are defined such that 
they satisfy the set of linear equations 


n-l 



k = 0 


i = 0, 1, ••• , n - 1 . 


( 20 ) 


The values of a t from Equation 19 are substituted into Equation 20, giving 



j =x 



0 ; 


i = 0, 1, ■ * • , n - J 


(21) 


For Equation 21 to be valid, it is sufficient for the x ; 's to be given by the roots of the polynomial, 


n “ 1 



k = 0 


(22) 


The bj *s can now be found using any n members of Equation 19. Thus, the solution of the non- 
linear problem is reduced to the solution of two linear systems plus the determination of the roots 
of an nth degree polynomial. From the roots of Equation 22, the values Uj , which divide the at- 
mosphere into slabs, are uniquely determined, using 


"J = TK) 1 ' < 23) 

The (AB). T s are obtained from the b - r s using the definition given in Equation 18. For the case of 
observations of I as a function of viewing angle the formulation is the same except that Equation 
17 is replaced by 


— - const + 'i ; 


l = 0, 1, ■ * * , 2n - 1 . 


King has applied the method to two model atmospheres with the Planck intensities given as 
analytic functions of optical depth, 


B(r) = 1 - e ’ T , 
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and 


B(t) = T e~ 2T 

Since the Laplace transforms of these functions are known, it was possible to compute values of 
the intensity- These intensity values were then used as "’observations” and the variable slab method 
applied to see if B(t) could be recovered. The results were encouraging. One somewhat obscure 
point in this method is the treatment of B(0), the Planck function corresponding to the temperature 
at the top of the atmosphere, which appears in the integration by parts form of the transfer Equa- 
tion 14. It has been pointed out by King (Reference 10) thatB(O) corresponds to the intensity seen 
at a viewing angle of 90° with a plane -parallel atmosphere, and he appears to treat it as a known 
quantity in his subsequent analysis. Since a viewing angle of 90° using a plane -parallel atmos- 
phere, has no physical validity it is not obvious how it can be handled in practice. 

The principal advantage of the method is the elimination of the need for arbitrarily choosing 
the divisions of the atmospheric slabs, gained at the expense of half of the observational data 
points available. The fact that the u } f s must be real and positive implies that the roots of Equa- 
tion 22 must be between zero and one in order to be physically acceptable, thus providing a poten- 
tial tool for error analysis. In order for the method to be of practical value, it is necessary for 
it to be generalized to include arbitrary band transmission. This problem will be considered in 
the next section. 

EXTENSION OF THE VARIABLE SLAB METHOD TO INCLUDE THE STRONG LINE 
APPROXIMATION FOR BAND TRANSMISSION 

Since nothing has yet been published concerning the extension of the variable slab method to 
arbitrary band transmission, a brief examination of the problem was made. It appears that in its 
present form the method requires a transmission with the functional form 

t(u, v) = e “*< v) * (u) (24) 

where and <£ are arbitrary functions. The exponential transmission of the strictly monochromatic 
case is, of course, of this form. The weak line approximation for both the statistical model and 
the Elsasser model is of a simple exponential form by Plass (Reference 14), thus satisfying Equa- 
tion 24. The strong line approximation for the statistical model is of the form 


t = e - (mu ) n 


(25) 


where in general m and n are functions of v. If the dependence of n on v is sufficiently small over 
the wave number interval being considered so that it can be replaced by an average value, then 
Equation 25 will be of the form of Equation 24. 
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In an effort to try the method on a more 
realistic atmospheric model and to compare 
it with another method, an attempt was made 
to fit the 15 m C0 2 transmission values, used 
by Wark (Reference 12) in his inversion cal- 
culations, to equations of the form of Equation 
25. Reasonable fits were obtained, being good 
toward the band center but becoming progres- 
sively less satisfactory away from the band 
center. Using four of Wark’s calculated inten- 
sities as ’'observational data," the variable 
slab method was applied, with the results shown 
in Figure 1. Wark's inversion calculation, 
which required only three data points, is shown 
for comparison. The horizontal lines show the 
slab divisions picked by Wark. 

In practice, the strong line approximation 
will probably not be satisfactory, especially for 

Figure 1 -Application of the variable slab method to a treating observations made in small wave num- 

model atmosphere. Calculated monochromatic intensities ber intervals. The Only entirely satisfactory 
at 690.7 cm -1 , 693.6 cm -1 , 697.2 cm” , and 705 cm -1 , j J J 

taken from Wark (Reference 12), were used as "observa- formulation would be one which permits the 

tions." The temperature profile used in the model is use of empirically given transmission functions, 

shown as the "sounding," and Work’s inversion calcula- 
tion is shown for comparison. 

APPLICATIONS TO THE MICROWAVE REGION 

While the above methods have been formulated with primarily the infrared spectral region in 
mind, the inversion problem can also be applied to the microwave region. Suggestions have been 
made for using the group of 0 2 absorption lines centered near 5 mm wavelength (References 15, 

16 and 17). Meeks has calculated weighting functions for this spectral region, using the Van Vleck- 
Weisskopf pressure-broadening theory in connection with a simplified model atmosphere. 

It would appear to be possible, at least in theory, to extract information on the thermal struc- 
ture at mesosphere altitudes, while in the infrared region the available information is confined 
primarily to the stratosphere. However, there are practical complications such as the necessity 
of including at high altitudes the effects of Zeeman splitting due to the earth’s magnetic field. 

THE QUESTION OF THE STABILITY OF SOLUTIONS 

Regardless of the method of solution used, the basic problem mathematically, as was pointed 
out above, is that of solving a Fredholm integral equation of the first kind 


f 


K(x, y) f(y) dy 


gO) 


(26) 
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where K(x, y) is the known kernel of the transform, f(y) is the unknown function and g(x) is a 
known function (in our case, the observational data). Now it can be shown that, for any kernel in- 
tegrable in the interval a to b, the following relation holds: 


lim f K( x, y) sin my dy - 0 . (27) 


Thus, the solution of Equation 26 tends to be unstable. Small uncertainties in g(x) can result in 
oscillatory solutions of large amplitude, instead of physically acceptable smooth functions. 

The right hand side of Equation 26 should actually be written g(x) + e(x) where e(x) is a small 
unknown error function. When written in this form, there is no unique solution for f(y). There- 
fore, it is necessary to supply conditions in addition to Equation 26 in order to pick out a unique, 
or at least a most probable solution from the infinite manifold of solutions available. One ap- 
proach is to set some conditions on the smoothness of the solution desired. This approach has 
been pursued by several investigators (References 6, 18 and 19). The main difficulty lies in the 
determination of how much smoothing should be introduced. Ideally one would like to employ just 
sufficient smoothing to eliminate the unwanted oscillations, while still retaining as much physical 
detail in the solution as possible. This would in effect represent the optimum solution obtainable 
from a given set of data. The inversion of Equation 26 is essentially a differentiation process, 
and as such it is much more sensitive to random errors than to systematic errors. 

It appears likely that whatever method of solution is adopted, it will probably be necessary to 
introduce some form of smoothing into the solution. 


CONCLUSIONS 

Of the various methods which have been proposed for solving the inversion problem, each 
seems to have its own peculiar advantages and disadvantages. The perturbation method has the 
advantage of utilizing a maximum amount of a priori knowledge of the atmospheric structure, at 
least when applied to the earths atmosphere. This method suffers from the disadvantage of re- 
quiring an arbitrary division of the atmosphere into slabs, as does the fixed slab method. The 
latter method, however, is perhaps best suited for direct smoothing techniques. Neither the vari- 
able slab method nor the analytic method require the undesirable arbitrary division of the atmos- 
phere. The latter method depends for its success on how well the thermal structure can be rep- 
resented by an expansion which is truncated after the first few terms. At the present time, the 
variable slab method suffers from a lack of complete generalization to include arbitrary band 
transmission. 

A considerable amount of work on the inversion problem remains to be done. Very little has 
been published in the way of error analysis and comparison of methods. Calculations should be 
made applying each method to the same realistic model atmospheres. The calculated intensities 
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serving as ’’observational data” can have errors artifically introduced into them so that the error 
propagation in each method can be studied. An effort should be made to determine whether the 
variable slab method can be satisfactorily extended to band transmission. The question of smooth- 
ing in relation to the variable slab method should be considered. 

Regardless of the method used, the problem is essentially that of solving a Fredholm integral 
equation of the first kind. Since the solutions of this problem are inherently unstable, it seems 
probable that some form of smoothing will have to be introduced into any method used. The amount 
of smoothing required will depend on the data and the magnitude of observational errors. The 
question then becomes one of whether the resulting smoothed solution retains sufficient physical 
information to be of practical value. The answer to this question will depend to some extent on 
the use to be made of the derived thermal structure, and it would seem to be profitable to examine 
the problem from this point of view also. 


(Manuscript received December 4, 1964) 
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Appendix A 


List of Symbols 

a. Constant defined by Equation 6 

a n Constant defined by Equation 12 

B(u ) Planck function averaged over band interval 
b j Constant defined by Equation 18 
b n Constant defined by Equation 12 
C. Constant defined by Equation 5 
C k Constant defined by Equation 20 
F. Orthogonal function 

f(y) Arbitrary function 

S(x) Arbitrary function 

H(u ) Chandrasekhar H- function 

1(0, fjL, v) Specific band intensity of radiation emerging from the top of the atmosphere in the 
direction () 

l(o, Observed outgoing radiation intensity 

I 0 (o, v.) Outgoing radiation intensity from model atmosphere 
K(x, y) Kernel of Fredholm integral equation of first kind 
k Mass absorption coefficient 

k a Roots of associated indicial equation 

l Constant defined by Equation 17 
m Constant defined by Equation 25 
n Constant defined by Equation 25 
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P Pressure 

S(r) Source function 

t Band transmission from level of normal absorber concentration u to top of the 

atmosphere 

u Normal absorber concentration (mass) in a vertical column 

x j Constant defined by Equation 18 

z Altitude 

a 2 Constant defined by Equation 18 

e(x) Small error function 

6 Zenith angle 

K Function equivalent to -at(u, v)/du 

M Cosine of zenith angle 

v Wave number 

ttF Radiant flux 

p Density of optically active gas 
r Normal optical depth in infrared 
4> Arbitrary function defined by Equation 24 

^ Arbitrary function defined by Equation 24 
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